;+
; NAME:
;       GRSMATCH
;
; PURPOSE:
;       Generate a spectrum for sources in the GRS overlap by
;       integrating over the source and subracting off a boundary spectrum.
;
; CATEGORY:
;       distance-omnibus

; CALLING SEQUENCE:
;       spectrum = GRSMATCH( struct [,V_STD=v_std][,SPECTRUM=spectrum]
;                           [,BGPS_ML=bgps_ml][,HI_DIRS=hi_dirs]
;                           [,BGPS_DIR=bgps_dir][,MAP_DIR=map_dir]
;                           [,OBJ_DIR=obj_dir][,/PLOT]
;                           [,ONSPEC=onspec][,BDRSPEC=bdrspec][,OBJMASK=objmask]
;                           [,BORDER=border][,X_ORIG=x_orig][,Y_ORIG=y_orig]
;                           [,X_BDR=x_bdr][,Y_BDR=y_bdr][,BHDR=bhdr][,R=r] )
;
; INPUTS:
;       STRUCT   -- BGPS source structure (see READ_BGPS_CSV), which
;                   includes source longitude and latitude
;       BGPS_ML  -- BGPS map locations IDL save file.  Defaults to:
;                   ./local/bgps_v102_map_locations.sav 
;
; OPTIONAL INPUTS:
;       GRS_DIR  -- Directory containing GRS spectral line cubes.
;       BGPS_DIR -- Directory containing the BGPS data
;       MAP_DIR  -- Directory containing the BGPS IMAGE files
;       OBJ_DIR  -- Directory containing the BGPS object files
;       R        -- Rind radius, in BGPS pixels, for averaging
;                   the "off spectrum".  [Default: 12]
;
; KEYWORD PARAMETERS:
;       PLOT     -- Plot the spectrum to the screen
;
; OUTPUTS:
;       SPECTRUM -- Integrated masked spectrum for the BGPS source
;
; OPTIONAL OUTPUTS:
;       V_STD    -- Velocity vector associated with the SPECTRUM
;       SPECTRUM -- Integrated masked spectrum for the BGPS source --
;                   may be different than RETURN'd spectrum if
;                   routine hits an early exit.
;       ONSPEC   -- "On-Source" HI spectrum weighted by BGPS emission.
;       BDRSPEC  -- "Border" HI spectrum, averaged over admissible
;                   border pixels
;       OBJMASK  -- BGPS catalog object mask -- pixels from BGPS
;                   assigned to this object.
;       BORDER   -- "Rind" around the BGPS source, excluding pixels
;                   assigned to other catalog objects.
;       X_ORIG   -- X-axis pixel numbers for the BGPS source
;       Y_ORIG   -- Y-axis pixel numbers for the BGPS source
;       X_BDR    -- X-axis pixel numbers for the border region
;       Y_BDR    -- Y-axis pixel numbers for the border region
;       BHDR     -- BGPS image header
;
; COMMON BLOCKS:
;       GRS_BLOCK
;
; MODIFICATION HISTORY:
;       Fri Jun 11 13:36:03 2010, erosolo <erosolo@>
;                 Documented
;       Thu Jul 22 11:14:30 2010, erosolo <erosolo@>
;		Updated to newest version.
;       Modified: 08/03/10, TPEB -- Extracted the spectrum creation
;                                   part of the original
;                                   prob_grsmatch.pro routine.
;       Modified: 06/15/11, TPEB -- Replaced DELVARX with UNDEFINE as the
;                                   means for clearing the data cubes
;                                   from memory before loading the next
;                                   one.  Fixed memory leak & improved
;                                   reliability.  Also, cleaned up
;                                   documentation.
;
;-

FUNCTION GRSMATCH, s, V_STD = v_std, SPECTRUM = spectrum, BGPS_ML = bgps_ml, $
                   GRS_DIR = grs_dir, BGPS_DIR = bgps_dir,$
                   OBJ_DIR = obj_dir, MAP_DIR = map_dir, PLOT=plot, $
                   onspec=onspec, bdrspec=bdrspec, objmask=objmask,$
                   border=border, x_orig=x_orig, y_orig=y_orig,$
                   x_bdr=x_bdr, y_bdr=y_bdr, bhdr=bhdr, R=r
  
  COMMON GRS_BLOCK
  
  ;; Set directory structure, if not specified in input
  if n_elements(grs_dir) eq 0 then grs_dir = './local/grs/'
  if n_elements(bgps_dir) eq 0 then bgps_dir = './local/bgps/'
  if n_elements(obj_dir) eq 0 then obj_dir = bgps_dir+'label/'
  if n_elements(map_dir) eq 0 then map_dir = bgps_dir+'maps/'
  
  
  ;; Define general parameters
  defsysv, '!MW', exists = exists
  if not exists then galactic_params 
  
  R0 = !mw.r0
  V0 = !mw.v0
  vs = !mw.vs    ;; Correction for HMSFR
  if n_elements(errvlsr) eq 0 then errvlsr = !mw.errvlsr
                                ; Assume a 7 km/s error in the
                                ; velocities

  ;; Define the velocity vector for use with the smoothed spectra
  v_std     = findgen(!MW.NVBINS)*!MW.DELTAV + !MW.VSTART
  spectrum  = fltarr(!MW.NVBINS)
  null_spec = spectrum
  onspec    = null_spec
  bdrspec   = null_spec
  
  
  ;; Read in Map Locations save file
  IF n_elements(bgps_ml) EQ 0 THEN $ 
     restore,file='./local/bgps_v102_map_locations.sav' ELSE $
        restore,file=bgps_ml
  
  
  hit = where(s.cnum eq bgps.cnum,ct)
  IF ct NE 1 THEN BEGIN
     message,'No match between distance structure and bgps SAVE file...',/inf
     RETURN, null_spec
  ENDIF
  bgps = bgps[hit]
  
  
  ;; Forward compatibility for ammonia spectral structures
  ;; flag = total(tag_names(bgps) eq 'BGPS_FILENAME')
  ;; if flag then filename = bgps.bgps_filename else $
  ;;    filename = bgps.filename
  ;; catroots = (filename)
  ;; slashpos = strpos(catroots, '/', /reverse_search)
  ;; catroots = strmid(catroots, slashpos[0]+1, 60)
  catroots = bgps.filename
  
  
  ;; Initialize structuring element
  IF n_elements(r) EQ 0 THEN r = 2
  
  elt = shift(dist(2*r+1,2*r+1),r,r) LE r
  
  ;; Check to see if this BGPS source lies within the GRS coverage
  IF bgps.glon LT 14 OR bgps.glon GT 60 THEN BEGIN
     message,'BGPS source lies outside GRS coverage...',/inf
     RETURN, null_spec
  ENDIF


  ;; Read in BGPS map and labelmap
  data = readfits(map_dir+catroots, bolo_hd, /SILENT)
  extast, bolo_hd, bolo_astrom
  sz = size(data)
  char = stregex(catroots, '_map')
  labelname = strmid(catroots, 0, char+1)+$
              'label'+strmid(catroots, char+1, 30)
  obj = readfits(obj_dir+labelname+'.gz', bhdr,/SILENT)
  
  
  ;; Determine where BOLOCAM points land in GRS data
  label = bgps.cnum
  objmask = obj eq label
  
  ind = where(objmask, ct)
  IF ct EQ 0 THEN BEGIN
     message,'No match in BGPS label map...',/inf
     return, null_spec
  ENDIF
  
  wts = data[ind]
  x_orig = ind mod sz[1]
  y_orig = ind / sz[1]
  
  border = (dilate(objmask, elt)-objmask)*(obj eq 0)
  ind2 = where(border,ctbdr)
  x_bdr = ind2 mod sz[1]
  y_bdr = ind2 / sz[1]
  
  xy2ad, x_orig, y_orig, bolo_astrom, glon, glat
  xy2ad, x_bdr, y_bdr, bolo_astrom, glonbdr, glatbdr
  
  wtorder = sort(glon)
  glon = glon[wtorder]
  glat = glat[wtorder]
  wts = wts[wtorder]
  lout = strcompress(string(round(glon)), /rem)
  filename = 'grs-'+lout+'-cube.fits'
  fn = grs_dir+filename
  
  wt2order = sort(glonbdr)
  glonbdr = glonbdr[wt2order]
  glatbdr = glatbdr[wt2order]
  lout2 = strcompress(string(round(glonbdr)), /rem)
  filename2 = 'grs-'+lout2+'-cube.fits'
  fn2 = grs_dir+filename2
  
  runspec = fltarr(!MW.NVBINS)
  bdrspec = fltarr(!MW.NVBINS)
  
  ;; Loop through all BGPS pixels assigned to this source
  FOR j = 0, n_elements(ind)-1 DO BEGIN
     IF ~ file_test(fn[j]) THEN CONTINUE
     
     ;; Pick GRS filename and load a new one if necessary
     IF fn[j] NE lastfn THEN BEGIN
        UNDEFINE,grs_data
        print,fn[j]
        grs_data = readfits(fn[j], hd)
     ENDIF
     lastfn = fn[j]
     extast, hd, astrom
     rdhd, hd, s = h
     sz_grs = size(grs_data)
     ad2xy, glon[j], glat[j], astrom, x, y

     ;; Look up spectrum
     if x lt 0 or x gt sz_grs[1]-1 or $
        y lt 0 or y gt sz_grs[2]-1 then begin
        wts[j] = 0
        continue
     endif
     spectrum = interpolate(grs_data,replicate(x,sz_grs[3]),$
                            replicate(y,sz_grs[3]),indgen(sz_grs[3]))
     
     ;; Interpolate to common velocity scale
     specout = interpol(spectrum, h.v, v_std)
     specout[where(v_std lt min(h.v) or v_std gt max(h.v))] = !values.f_nan
     runspec = runspec+specout*wts[j]

  ENDFOR   ;; End of loop for ON-SOURCE spectra extraction
  
  
  ;; Loop through all BORDER pixels
  FOR j = 0,n_elements(ind2)-1 DO BEGIN
     IF ~ file_test(fn2[j]) THEN CONTINUE
     
     ;; Pick GRS filename and load a new one if necessary
     IF fn2[j] NE lastfn THEN BEGIN
        UNDEFINE,grs_data
        print,fn2[j]
        grs_data = readfits(fn2[j], hd)
     ENDIF
     lastfn = fn2[j]
     extast, hd, astrom
     rdhd, hd, s = h
     sz_grs = size(grs_data)
     ad2xy, glonbdr[j], glatbdr[j], astrom, x, y
     
     ;; Look up spectrum
     if x lt 0 or x gt sz_grs[1]-1 or $
        y lt 0 or y gt sz_grs[2]-1 then begin
        ctbdr = ctbdr - 1
        continue
     endif
     spectrum = interpolate(grs_data,replicate(x,sz_grs[3]),$
                            replicate(y,sz_grs[3]),indgen(sz_grs[3]))
     
     ;; Interpolate to common velocity scale
     specout = interpol(spectrum, h.v, v_std)
     specout[where(v_std lt min(h.v) or v_std gt max(h.v))] = !values.f_nan
     bdrspec = bdrspec+specout
     
  ENDFOR   ;; End of loop for OFF-SOURCE spectra extraction
  
  ;; Filter the spectra
  filter = savgol(17,17,0,6)
  
  specout = runspec/total(wts)
  spectrum = convol(specout,filter,/edge_trun,/nan) >0
  
  bdrspec = bdrspec/ctbdr
  bdrspec = convol(bdrspec,filter,/edge_trun,/nan) >0

  
  ;; Create the Final Spectrum     
  onspec = spectrum
  emmask = spectrum GT 4*mad(spectrum)
  spectrum = ((spectrum - bdrspec)>0)*emmask
  
  
  ;; Plot, if desired
  IF KEYWORD_SET( plot ) THEN BEGIN
     window,1
     plot,v_std,onspec,/xst,xtitle='V!DLSR!N (km/s)',$
          ytitle='!6T!DA!N (K)',linestyle=2
     oplot,v_std,bdrspec,ps=10,color='0000ff'x,linestyle=2
     oplot,v_std,spectrum,ps=10,color='ffff00'x

     v_tan = V0*(1. - sin(s.glon_peak * !dtor))
     plots,v_tan,!y.crange[0]
     plots,v_tan,!y.crange[1],/continue,linestyle=2
     print,!x.crange
     print,n_e(v_std),min(v_std),max(v_std)

  ENDIF
  
  ;; Square the spectrum  -- WHY?
  spectrum = spectrum^2

  IF total(spectrum) EQ 0 THEN RETURN,null_spec
  
  RETURN,spectrum
  
END

